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Abstract: 

We construct a new global potential energy hypersurface for H 2 + H 2 and 
perform quasiclassical trajectory calculations using the resonance complex theory 
and energy transfer mechanism to estimate the rate of three body recombination 
over the temperature range 100 - 5000 K. The new potential is a faithful rep- 
resentation of ab initio electronic structure calculations, is unchanged under the 
operation of exchanging H atoms, and reproduces the accurate Hz potential as one 
H atom is pulled away. Included in the fitting procedure are geometries expected 
to be important when one H 2 is near or above the dissociation limit. The dynamics 
calculations explicitly include the motion of all four atoms and are performed effi- 
ciently using a vectorized variable-stepsize integrator. The predicted rate constants 
are approximately a factor of two smaller than experimental estimates over a broad 
temperature range. 


Mailing address: NASA Ames Research Center, Moffett Field, CA, 94035 


I. Introduction 

We have begun a detailed study of the calculation of recombination rate con- 
stants for hydrogen containing compounds in the gas phase. This is motivated by 
the need for these rate constants as input into modeling studies of hydrogen com- 
bustion processes. These processes in turn are important for the elucidation and 
design of a wide variety of hydrogen containing systems. 

The initial system being considered is the recombination of hydrogen atoms 
to form H 2 in the presence of excess thermal if 2 . This system has the advantage 
that it is one of the simplest possible and also there exist several experimental 
measurements which can be used to monitor the reliability of the theoretical pre- 
dictions. The available experimental studies have been reviewed [l] and values for 
the recombination rate constants have been recommended over the temperature 
range 50-5000 K. The availability of experimental data will be important because 
the theoretical models used are not completely developed and the ultimate goal is 
to reliably predict rates for other systems or under conditions where experimental 
results are not available. 

There have been several previous theoretical studies of gas phase recombination 
processes and different models and methods have been developed[2-6]. The particu- 
lar model used in the present study was formulated by Roberts et al. [4] whereby the 
recombination rate is calculated from the rate of stabilization of quasibound reso- 
nant states by collisions with a third body. Of the available models, this resonance 
complex theory is the only one which takes into account the quantum mechanical 
nature of the metastable states involved in the recombination process. This model 
has been utilized in conjunction with the quasiclassical trajectory method to calcu- 
late recombination rates for hydrogen using several third bodies( if, if 2 , He, and 
Ar in Ref. [7] and H in Ref.[8]). The present study closely follows this previous work 
but also includes several extensions and improvements. 

The computational steps involved in the resonance complex theory are twofold. 
First of all, it is necessary to identify and characterize the important resonant states 
and secondly the rates of stabilization of these resonant states must be determined. 
For the present work we consider only the energy transfer mechanism in which 
case the resonant states are the metastable states of if 2 - Previous workers have 
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well characterized these metastable states[9-ll], thus in the present report we axe 
concerned with the four body problem of collisions of thermal Hi with metastable 
Hi. 

The previous calculations for Hi as a third body [7] made several approxima- 
tions. Probably the most severe limitation was the use of empirical potential energy 
functions. Another possible source of error was the assumption that the thermal 
Hi could be treated as a structureless particle. In the present work we eliminate 
these criticisms by constructing a new potential energy function which is based on 
ab initio electronic structure calculations and accurately treat the full four body 
dynamics. In addition, we use more accurate metastable state information and 
extend the results to higher temperatures. 

The new potential is described and discussed in Sect. II, the dynamical meth- 
ods are described in Sect. Ill, and our results are presented and discussed in Sect. 
IV. Finally in Sect. V we present the conclusions of this study. 

All our calculations are carried out using hartree atomic units where the unit 
of energy is denoted Eh and 1 Eh — 219474.7 cm -1 = 627.5095 kcal/mol = 4.3598 
xlO~ 18 J, the unit of length is the bohr and is denoted a 0 and 1 a 0 = 0.5291771 
xlO -10 m, the mass unit is the electron mass which is denoted m e and 1 m e = 
5.485804 XlO -4 a.m.u., and the time unit is equal to 2.418884 xl0~ 17 seconds. 

II. Potentials 

First consider the Hi potential. The quasibound states of H 2 were studied 
for several potentials in Ref. [ll] and because of the variation of the results seen 
there, it is desirable to employ the most accurate potential curve possible. The 
most accurate resonance parameters determined in Ref. [ll] used a potential which 
included estimates of nonadiabatic effects. However this potential is not suitable for 
classical dynamics because it can not be represented as a single potential curve, thus 
instead we use the potential which includes all corrections except the nonadiabatic 
correction. The potential we will use is a representation of the one called Ad in Ref. 
[ll]. This potential is based on the most accurate Born-Oppenheimer potential aug- 
mented by relativistic, radiative and nuclear motion corrections. The representation 
of the Ad potential used in Ref. [11] was a combination of a piecewise fifth order 
polynomial fit to the Born-Oppenheimer potential and its first derivative and cubic 
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spline representations of the various corrections to the Born-Oppenheimer potential. 
In addition, extrapolations using various functional forms for small and large bond 
length are involved when evaluating that potential. For the present calculations, it 
is desirable that the potential be represented in a form which avoids decisions in 
order to take full advantage of the vector processing capabilities available on class 
VI supercomputers. Thus we search for an alternate analytic representation. 

We begin by writing the new ground state potential V H H in the form 

V H H = V SR {r)+V LR (r ), ( 1 ) 

where r is the bond length, V SR is a short ranged contribution and V LR is the long 
range form of the potential, which is represented as 

^ iR = -C 6 /(r 6 + rS)-C s /(r 4 + r3) J -C l0 /(r J + r2) 5 . (2) 

We take the dispersion parameters Cq, Cs, and C\o from previous calculations(l2,13] 
and treat ro as an adjustable parameter. After some experimentation, we chose to 
represent V SR as 

= D{ezp\—a(r)] - l} 2 - D, (3) 

i.e. a Morse curve with a r dependent range parameter. Taking cognizance of the 
singularity at r = 0, we write 

N 

a = r _a ^a t r t . (4) 

t=0 

The fitting procedure consists of first guessing values for ro, s and N, then using Eqs. 
(1) and (2) to determine V SR . The parameter D is then fixed by approximating a 
by a(r - ri) and determining D, a and ri so that Vsr is reproduced at the three 
values of r closest to the equilibrium separation of if 2 . Once D is known, Eq. (3) 
can be inverted to give a at the input data points. The remaining parameters are 
then determined by linear least squares. This procedure is repeated as a function 
of ro in order to minimize the rms error produced in the fit to a. The data used in 
the fit are the values of the Ad potential from Ref. [11] evaluated at the distances 
where the Born-Oppenheimer potential was calculated[l4|. 


4 



An important problem for least squares fitting procedures is the method of 
weighting the various data points. For example, in the above procedure, an equally 
weighted least squares fit would tend to fit most accurately the values of a for large 
and small r since the magnitude of a is largest for these extremes. However, it is 
probably desirable to fit a most accurately for intermediate r, since that is where 
the wavefunction is likely to be largest. We accomplish this goal by using a weight 
function which is peaked in the vicinity of the equilibrium separation. The particular 
function used is the square of the ground state Morse oscillator wavefunction for the 
Morse parameters jD, a, and rj determined above. In addition to these parameters 
a mass is required in the weight function - this parameter determines the width of 
the weight function. Empirically a value of 10 m e was determined to be suitable. 

In Fig. 1 the curves a and r*a for s=l are given. Both curves are very smooth 
and appear to be much easier to fit than the potential itself. The curve for ra looks 
approximately quadratic and so it might be expected that a low order fit would 
probably be quite adequate. However, because the potential for H 2 is known very 
accurately, it is necessary to use a high value of N to obtain satisfactory results. 
The final fit is obtained using s=l and iV=16 in Eq. (4). It should be noted that an 
important reason for the approximate linearity of a for large r is the separation of 
the long range and short range form of the potential. It would be interesting to see 
if these techniques would be useful for representing the potential of other diatomics. 

An important problem when using high order polynomials is to make sure 
that they do not do strange things for bond lengths not included in the fit. This 
caused many trial fits to be rejected. The current fit to ra has only one zero for 
r real and r > 0 (the zero is near the equilibrium separation) and also d\ra}/dr 
has only one zero for r real and r > 0. The overall rms error of the fit is 2.78 
mEh, however the error decreases rapidly as r increases and the rms error for all 
points with r > 1 a 0 is only 0.52 [iEh- We have also computed the energy levels 
of all bound and quasibound states of H 2 using the WKB method for the current 
analytic representation and the original Ad potential and they agree very well. The 
rms difference is only 0.16 cm -1 while the maximum difference is only 0.30 cm -1 . 
Thus the current analytic representation satisfies our computational demands of no 
decisions for its evaluation without sacrificing accuracy. The parameters for this 
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representation are given in Table I. It should be noted that because the potential 
includes corrections to the Bom-Oppenheimer potential which are mass dependent, 
it is only valid for the isotope of H with mass 1 a.m.u. 

We now turn to the full H 4 potential. There have been several potentials 
proposed for the H 4 system, however few are applicable when one or more hydrogens 
are not at their equilibrium displacement and of these, fewer still are of documented 
accuracy when a hydrogen is stretched far from equilibrium. However, it is just these 
geometries which are expected to be important for the present application. Thus we 
will construct a new potential which concentrates on this aspect. The data we use to 
produce the new potential is the result of ab initio electronic structure calculations. 
There have been a large number of studies of the electronic structure of H 4 , again 
mostly for hydrogen at its equilibrium separation. Because the classical equations 
of motion use as input only the gradient of the potential and not the potential 
itself, it will be important to accurately know the dependent of the potential on 
all degrees of freedom including bound length. Thus we perform new ab initio 
electronic structure calculations with these requirements in mind. 

The electronic structure calculations where carried out using the MOLECULE- 
SWEDEN[15] codes with a gaussian basis set. The basis set consisted of 8 s func- 
tions contracted to 4 functions, 2 p functions and 1 d function per atom. The s 
function exponential parameters where taken from Ref. [16] and the contraction 
coefficients where determined from a calculation on the hydrogen atom. The p and 
d function exponential parameters where taken from Ref. [17]. The s components 
of the d functions were not included in the calculation giving a grand total of 60 
contracted basis functions. With this atomic orbital basis, molecular orbitals where 
determined by performing a CASSCF calculation which included 4 orbitals in the 
active space. The energy produced by this calculation is denoted E CAS . Using this 
molecular orbital basis, a multireference configuration interaction calculation was 
carried out including all single and double excitations out of the possible CAS ref- 
erence configurations. The energy produced by this calculation is denoted E MRCI . 
All calculations include all four H atoms to minimize problems of size consistency. 
With this basis set and correlation treatment, the difference in E MRCI between 
two i ?2 molecules with bond length 1.401 a 0 separated by 20 a 0 and four H atoms 
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no closer than 20 a 0 to each other is —0.345698 Eh , which is 99.1 % of twice the 
accurate value of the Born-Oppenheimer D t . 

Calculations were carried out for a total of 92 different geometries, all using 
C* 2 u symmetry, and the results are given in Tables II- VI. The number of geometries 
may seem low, but due to the high symmetry of this system and the method used 
to represent the potential energy hypersurface, these points should be sufficient for 
our purposes. 

The ab initio electronic structure calculations we have carried out are expected 
to be accurate except in the vicinity of the Van der Waals minimum where much 
larger basis sets and more accurate correlation treatments are required(l8j. This 
should not be a problem for two reasons. First of all it is possible to represent the 
potential in the vicinity of the Van der Waals minimum by using data describing the 
asymptotic form of the interaction potential, and secondly, the energies sampled by 
our dynamics calculations will be larger than the Van der Waals well depth, thus 
small to moderate errors there should not have a large effect on the calculated rate 
constants. 

To represent the long range potential not given accurately by our electronic 
structure calculations, we will include a function in our analytic representation 
which has the proper asymptotic form. Most of the parameters of this function are 
determined not by the present electronic structure calculations but rather by other, 
more accurate, calculations. In particular the long range potential is represented as 

Va£+CD = ^2 V q 1 % f i{ r l^ r 2,R)yq l q 2 ^i^2,<P), ( 5 ) 

q iq^n 

where AB A CD denotes a particular way of drawing the bonds between four indis- 
tinguishable hydrogen atoms, r i , r 2 , R, 6\, 62 , fa, and fa are the jacobi coordinates 
for AB -tCD, and <p = 4>\ — fa. The jacobi coordinates are defined so that r, is the 
length of the vector r t connecting the atoms in molecule i, R is the length of the 
vector R connecting the centers of mass of the two molecules, cos0 t - = r, • R/Rn, 
and cos<j> = (r\ x R) • (r 2 x R)/\r\ x R \\?2 x R\. The angular functions are given 
by [19] r " r '" 

4tt 

= [ 2(1 + 6 o)] 1 / 2 

( 6 ) 
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where is a spherical harmonic. Because of homonuclear symmetry, only even 
values of the indicies 91 , and 92 appear in Eq. (5). In Eq.(5), we will only include 
the terms with q\ + 92 < 2. The terms with 91+92 <2 include dispersion and 
induction contributions and the terms with 9i + 92 = 2 also include the quadrapole- 
quadrapole interaction. We will include the dominant contributions for these terms 
and write 

V™, = \ 91 + 92 <2; (7) 

l j2’'”'<?(ri)<?(r 2 )/[J? 5 +<i is (tf) 5 ], 71+92=2, 

where C| l92M and Cg 1,2M are dispersion coefficients, d LR is a damping function, Q 
is the quadrapole moment, and 


f 3/10, p = 0; 

Q qiq3tl = < V2/5, n — l; 

[y/2/20, fj, — 2. 


( 8 ) 


Everything except the damping function is obtainable from calculations by other 
workers. The damping function is determined by smoothly blending V^b+cd 
the rest of the potential to approximately reproduce the Van der Waals minimum. 
In particular we use 

d LR = df R exp(-d% R R). (9) 

For fitting purposes it is convenient to independently vary the quantities d% R and 
d LR (R = 6.5 a 0 ) rather than d% R and d{ R . Then the parameter d\ R primarily 
controls the size of the long range potential for geometries which are mostly repulsive 
and d LR (R = 6.5 a 0 ) controls the depth of the Van der Waals well. 

The dependence of the quadrapole moment on bond length is parameterized 
by fitting the calculated values of Q from Ref. [20] over the range r = 0.4 to 4a 0 . 
We choose to represent it as 


Q = r 2 (a Q + b Q r + c Q r 2 )exp{-{d Q r) 2 ], cr. - ( 10 ) 

with the parameters determined by equally weighted least squares. The parameters 
are given in Table VII. The rms error for this fit is 4.2xl0~ 3 o.u. 
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The dependence of the dispersion parameters on bond length is not as well 
known, so it is necessary to resort to a more approximate procedure to parameterize 
this dependence. We do this by combining the values of the parameters for r x = 
f 2 = 1.449a 0 quoted in Ref. [21] along with the dispersion coefficients for He + H 2 
from Ref. [18], which are given for three bond lengths. In particular, we write 

= [(291 + 1)(2, 2 + l)]- ,/a CS‘(ri)C«>(r,)/C„(ff«), (11) 

where C£ is the dispersion coefficient for He 4- H 2 multiplying the Legendre poly- 
nomial of order q, and C n {He) is an effective dispersion coefficient for He + He 
interactions, chosen so that equality holds in Eq. (11) when r x = r 2 = 1.449a,,. 
The square root factor in Eq. (11) arises because of the normalization factors of 
the spherical harmonics in the angular function l/ 9l<j2M not present in an Legendre 
polynomial expansion. The bond length dependence of the Cl is parameterized by 
writing 

Cl = r m "(a* + b q n r)exp[-{d Q r) 2 }, (12) 

i.e. we assume that the r dependence of C% is similar to that of the quadrapole 
moment. The parameters for this equation are determined by equally weighted least 
squares. In Table VII we give the parameters for the functions Cl which are the 
Cl rescaled so that the equation 

Cl^ = Cl'{r x )Cl'{r2) (13) 


is satisfied. 

Equation (11) is an approximation, and it is interesting to estimate its accuracy 
by comparing the C n (He) to accurate values. For the case n = 6, we obtain values of 
Cq(H e) which range from 1.24 to 1.33 a.u., while the accurate value is 1.47±0.1 a.u. 
[22], and for n = 8, 14.9 to 16.8 a.u., as compared to 14.0±0.2 a.u. [22] Thus we 
see that Eq. (11) holds to about 20 %, which should be accurate enough for the 
present calculations. 

We now turn to the representation of the repulsive part of the potential. One 
of the difficulties of producing an analytic representation for all atoms near together 
is the high symmetry in the H 4 system. This is especially true if the representation 
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needs to behave properly when H atoms are exchanged. When both molecules are 
near their equilibrium separation, it is convenient to expand the potential in terms 
of the jacobi coordinates r l5 r 2 , R, 6 i, 0 2 , and <i>, as was done for the long range 
potential. However when R and the r,- become similar in size, it is not always 
clear how to drawn in the bonds, i.e. it may be ambiguous whether the system 
is best described by AB + CD, AC + 1 3D or AD + BC, where A-D are formally 
distinguishable H atoms. Thus a given set of jacobi coordinates are not suitable for 
a global representation of the potential. In the present situation we will handle this 
problem by first determining a representation for small vibrational displacements 
where the identity of the bonding and nonbonding pairs is clear, then switch to a 
more suitable representation for other configurations. 

We begin our representation by extending the techniques of Ref. [23]. We first 
define the interaction potential to be the difference in energy between a geometry 
having particular values of r l5 r 2 , R, 9\, 0 2 , <t> and a geometry with the same values 
of r i and r 2 but with R = 20 o 0 . The total potential is the sum of the interaction 
potential and the bonding pair potentials, V HH . The interaction potentials from 
the electronic structure calculations use E MRCI and are denoted V MRCI . 

The first ingredient in the analytic representation is a nonbonding pairwise 
potential called Vp. This potential is determined by considering the interaction 
potential of the geometries having 0\ = 0 2 = cf> = 7r/2. This is called the crossed 
geometry and its energies are given in Table VI. For these geometries, when both 
molecules have the same bond length, all of the nonbonding pairs are separated 
by equal distances, thus the nonbonding pairwise potential can be defined as one 
quarter of the interaction potential. To make this determination unique, only the 
eight geometries in Table VI having bond lengths equal to 1.401 a 0 are used. This 
determines Vp only for eight distances, and we extrapolate and interpolate V P by 
fitting it to the functional form 


Vp = Aptxp[— Bpx ( ~' F ), (14) 

with Ap, Bp, and Cp parameters. This fit is performed using nonlinear least 
squares with each point weighted by the amount (\V MRCI /4\ + lOmi^) -1 , where 
yMRCi j s the interaction potential calculated from Table VI for that distance. Thus 
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the approximate relative error is minimized. The fit is quite good, giving as the 
minimized error 3.7 x 10 -2 and an absolute rms error of 10 nE\- The parameters for 
this fit are given in Table VII. It should be noted that Vp is positive everywhere. In 
Ref. [23], the nonbonding pairwise potential was generated at arbitrary distances by 
cubic spline interpolation, however we use Eq. (14) instead because of the smaller 
number of input data points and for the vectorization considerations mentioned 
above. 

The nonbonding pairwise potential provides a reasonable zeroth order repre- 
sentation of the interaction potential, but it is not sufficiently accurate, so we make 
corrections to it. For small vibrational displacements, we write 


V AB+CD ={Vp[R A c) + Vp(Rxd) + Vp(Rbc) + Vp{RBD)]F c {rur2,R,6i,02,<t>) 


+ ^AB+CD( r l) r 2)-JMl>02)<£) 


(15) 


where is the interaction potential for AB and CD close to their equilibrium 

displacements, Rij is the distance between atom i and atom j , F c is a multiplicative 
correction function, and V%jj +CD is the long range potential defined above. In Ref. 
[23], a small factor was included to shift the potentials to ensure that F c was always 
positive and always multiplied a positive function, however that is not required here 
because in the present case the long range attractive part of the potential is treated 
separately so that Vp and the differences between interaction potentials and the 
long range potentials included in the fitting procedure are automatically positive. 

We now expand the function F c as 


F‘ = £ f^Jr u r 2 ,R)y , *), (16) 

<7i<72M 

where the coefficients are to be determined and the angular functions are 

given by Eq. (6) above. To determine the , we use the difference between the 

interaction potential and the long range potential for r, = 1.0 or 1.401a o , R = 3.0, 
3.5, 4.0, 4.5, or 5.0 a 0 , and [ 81 , 62 , <f>] = [0,0,0] (linear), [tt/2,0,0] (T shaped), 
[0,7t/2,0] (T shaped), [tt/2, tt/ 2,0] (parallel) or [tt/2, tt/2, 7t/2] (crossed), a total of 
100 energies. This plus Eq. (15) gives F c . The terms {qi,q 2 ,fj)= (0,0,0), (2,0,0), 
(0,2,0), (2,2,0), and (2,2,2) are used in Eq. (16), and this equation is inverted to 
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give the The term (91, 92> m )=(2,2,1) is not included in the fitting procedure 

because it is zero for all the geometries for which we have calculated ab initio data. 
To generate the coefficients f$ iqalt for arbitrary values of r» and R, they are fit by 
linear least squares to the form 


1 1 


f c = V' Y' V' c iik 

t=o y=o Jt=o 


r'r'fl 4 . 


( 17 ) 


All points are equally weighted in this fit and typically rms errors of about 0.3 mEh 
are obtained. The final parameters are given in Table VIII. 

It should be noted that the restriction of the angular functions in Eq. (16) 
to + 92 < 2 does not imply a similar restriction on the potential VJjj+ci? - the 
presence of the pairwise potential V P ensures that higher order anisotropies exist. 

The analytic representation of Eq. (15) works quite well when the bond lengths 
are not too large, but becomes less accurate for large displacements from equilib- 
rium. For a typical example, see Fig. (2). We now define a large vibrational 
displacement potential which is more accurate for large bond lengths. A potential 
valid for both small and large displacements is obtained by smoothly switching from 
one to the other. 

A subset of geometries where one bond length is very large is the situation 
where the system is most naturally described as Hz + H , thus the large vibrational 
displacement potential will be forced to satisfy this limit. In particular, when one 
H is pulled away, the large vibrational displacement potential will reproduce a 
modification of the accurate LSTH Hz potential[24]. The LSTH total potential is 
written in the form 


Y LSTH 


V Lon 


+ V" SB 


' HHHi 


(18) 


where Vfi° n is obtained from the london equation: 


y«r = E 2,- - d Ei-* - -y 2 > 1/2 . 

*=1 t >j 

where 

£.' = 5[‘*W + 3 £(*<)!. 


(19) 


( 20 ) 
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( 21 ) 


J< = - 3 W 

Ri is one of the three interpair distances, 1 E is the H 2 singlet potential, 3 E is an 
effective H 2 triplet potential, and V^hh a correction function for nonsymmetric 
and bent geometries. Both V/j(° n and V^hh are symmetric with respect to the 
interchange of any two H atoms. The simplest way to extend V LSTH to H A is to 
write 

6 , 
i=l 

( 22 ) 

where fit and Ji are defined by Eqs. (20) and (21) with interpair distances R\ — 
Rab , - Rac, R3 = Rad, Ra = Rbc, Rs = Rbd, and Rq = Rcd, and to write 


(23) 


and 

v«. = v£" + Vgf B . (24) 

In practice, this is modified in three ways. First of all, we replace the singlet function 
l E of Eqs. (20) and (21) with the H 2 potential Vjjh. Secondly we slightly modify 
the effective triplet potential 3 E so that it is always greater than Vjjh by writing 

3 £= 3 s{R)V£>f{R)+ 3 E{R), (25) 


where 

3 s = exp{— 3 ai? 2 ], (26) 

and 

Vffi = B{exp[-a(fl)] + l} 2 - D + V lr (R), (27) 

D, a, and Vlr defined at the beginning of this section. The parameter 3 a is 
given in Table VII and is determined to approximately minimize the effect of this 
change while ensuring the inequality 3 E > V HH . These two changes do not signifi- 
cantly alter the LSTH potential for geometries probed by most collisions. The final 
modification is to scale V^ SB in Eq. (24) by a factor which goes to unity as the 
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configuration goes to Hz + H. Thus the large vibration displacement total potential 
is 

V LD = Vk 0 A n + S{p)Vgf B , (28) 

where V£° n is given by Eq. (22) with the replacements V HH and Z E for l E and 
3 E in Eqs. (20) and (21), and the scaling function is given by 

S = 1/{1 + exp[-0(p - />„)]}, (29) 

where 

(30) 

t=l 

The quantities 0 and p Q are parameters. The V LD potential is symmetric with 
respect to the interchange of any two H atoms. The motivation of including the 
factor S can be seen in Fig. (2) where the interaction potentials derived from V/f° n , 
and V LD are shown. In comparison to the ab initio data, the Vjj° n potential 
is much too small and the potential is much too large. 

The small and large vibrational displacement potentials are now combined by 
defining the quantity 

Vab+cd = ^ D (ri+T2m^i'i D +V HH (r l )+V HH (r 2 )\+[l-3 D (r l +r 1 )}V LD , (31) 
where the switch s D is given by 

s D = ll{l + exp[a D {r 1 +r 2 -r D )\}, (32) 

where a D and r D are parameters. 

The final linear parameters in Eq. (17) are determined by approximately op- 
timizing the nonlinear parameters 0 , p 0 , a D , d% R , and d LR (R = 6.5 o 0 ). The 
parameter r D is fixed at a value of 5.4 a c , which causes s D to equal 1/2 when r 2 
equals 4 a 0 [see Fig. (2)]. The three nonlinear parameters 0 , p 0 and a D are adjusted 
to approximately minimize the weighted error in the interaction potential for all ge- 
ometries having r\ or r 2 greater than 2 a 0 . The weight function used when rj or 
r 2 > 2a 0 is (\V MRCI \ + 1 mEh) -1 , where V MRCI is the ab initio interaction poten- 
tial at these geometries. The parameters d% R and d LR [R = 6.5 o 0 ) are adjusted to 


14 


minimize the error in the points with r* small and to produce good agreement with 
the potential expansion coefficients in the vicinity of the Van der Waals minimum 
given in Ref. [25] . The final parameters are given in Tables VII and VIII. The rms 
error for the fit to VXb+cd 0-29 mEh and the weighted rms error for the points 
with r\ or r 2 greater than 2 a 0 is 1.8xl0~ 3 . 

The quality of the Van der Waals well can be estimated from Fig. (3). Here 
the first two functions \{qi,q 2 ,n) = (0,0,0), (2,0,0)] of an expansion like Eq. (5) 
of Vab+cd is shown along with the data used to adjust the parameter d LR (R = 
6.5 a 0 ). The agreement is quite good. To be consistent with any systematic errors 
introduced in calculating the spherical average of the potential, we used the same 
method as in Ref. [25], namely that the interaction potential at the four geometries 
(linear, T-shaped, parallel, and crossed) and the difference between the interaction 
potential for a trapezoidal and parallelogram geometry is fit to a five term angular 
expansion. In addition, this figure shows the results of accurate calculations of the 
expansion functions. These were obtained by using converged numerical quadrature 
and the orthogonality of the l/ 9l( j 3M . The two methods of calculating the expansion 
functions do not differ greatly. Finally, this figure show the accurate expansion 
functions for the potential denoted DHR by Ref. [7]. These will be compared in 
detail in Sect. IV. 

It should be noted that the angular expansion we use is different than that 
used in Ref. [25], however the two expansions are equivalent and are related by a 
nonsingular transformation. 

In addition to the expansion functions, we have computed the geometry of 
the Van der Waals minimum and it has R= 6.27, r\ — r 2 = 1.402 a 0 , 6\ = ?r/2, 
0 2 = <t> = 0, and has a dissociation energy of 0.14 kcal/mol. 

The remaining facet of the fit concerns the drawing of the bonds, i.e. whether 
the system is best described as AB + CD, AC + BD , or AD + BC. This ambiguity 
is resolved by defining a weight function for each of the three ways of drawing the 
bonds which picks out the appropriate way. To do this, define the quantity r) A B+cD 
as 

r}AB+CD = [(rf 2 + r 2 2 ) l/2 -R]"\ (33) 

where r 1? r 2 and R are a subset of the jacobi coordinates for drawing the bonds 
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as AB + CD. Of the three ways of drawing the bonds, the minimum value of 
Vab+cDi Vac+bd and tjad+bc will likely indicate the most physical situation. 
The motivation for the functional form of Eq. (33) is as follows: in general the 
bonds should be drawn to minimize the individual bond lengths and maximize R. 
A differentiable expression for mm(ri, rj) is obtained by using [r^ n + r^ n ] -1//n and 
taking the limit as n goes to oo, however the simplest case of n = 2 seems adequate 
for the present application, thus tjab+cd Is approximately equal to min{r\ y ri) / R. 
The motivation for maximizing R can be seen by considering a linear arrangement 
of the atoms in the order ABCD. If the atoms B and C were closer together than 
A and B or C and D , then a rule based only on mm(ri, r 2 ) would select AD + BC 
as the appropriate way to draw the bonds. However if the distance between A and 
B and the distance between C and D were equal, R for the arrangement AD + BC 
would be zero. Thus maximizing R also would, for this example, select AB + CD 
as the appropriate way to draw the bonds. Combining these ideas, we take the final 
potential to be 

V =[S a (t]a.b+cd)Vab+cd + S a (tjac+bd)^ac+bd + S A {tIad+bc)Vad+bc]/ 
\S A {rfAB+CD ) + S a {ti ac +bd) + S a (tiad+bc)], 

(34) 

where 

S A (rj) = exp[-P A ri], (35) 

/ 3 A a parameter we take to be 20. This final potential is symmetric with respect to 
interchanging any two H atoms. 

An interesting quantity for a potential such as the current one which allows the 
possibility of hydrogen atom exchange, is what is the barrier for such an exchange. 
If the transition state for exchange is the square planer geometry, then the current 
potential predicts a barrier of 144 kcal/mol with a side of length 2.7 a 0 . This 
compares favorably with the ab initio value of 140 kcal/mol at 2.4 a 0 [26,27]. 

One unfortunate aspect of the present potential is that it is fairly time con- 
suming to evaluate. In the trajectory calculations we require the gradient of the 
potential, and we calculate the gradient with respect to the six atom-atom distances 
and transform those derivatives using the chain rule to the form required by the 
equations of motion discussed below. The majority of the time spent evaluating 
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the gradient is spent determining the factor d^/dRij for the three different sets of 
jacobi coordinates. These factors are required for any potential which is symmetric 
with respect to exchanging any two hydrogen atoms and has the correct long range 
quadrapole-quadrapole interaction. Thus there is nothing special about the present 
potential which makes it expensive to evaluate. 

III. Dynamics Calculations 

We follow Roberts et al. [4] and calculate the rate constant k 3B defined by 

d{H 2 )/dt = k ZB \H} 2 \M), (36) 

for the process 

H + H + M->H 2 +M , (37) 

by modeling the recombination process as occurring via the steps 

H + H~H;{i), (38) 

and 

H* 2 {i)+M - H 2 {n) + M, (39) 

where M is thermal H 2 , H 2 is a quasibound state of H 2 , i labels a particular 
quasibound state and n labels a particular bound state. This is the energy transfer 
mechanism. If the rate determining step is the process of Eq. (39) and the process 
of Eq. (38) is at equilibrium, then the rate constant can be written 

* 3B = E K i,m k Lm = E k,w)H cn. ( 4 °) 

i } n i 

where K\ q is the equilibrium constant for Eq. (38), k{ n is the rate constant for Eq. 
(39), k{ is k{ n summed over n, and T is the temperature. 

If Eq. (38) for a particular quasibound state is not at equilibrium, then k 3B as 
calculated via Eq. (40) would be too large if that quasibound state was included in 
the sum. Thus we estimate nonequilibrium effects by simply neglecting all contri- 
butions to the overall rate due to that quasibound state[4]. Hence the quasibound 
states are partitioned into two groups, those that are approximately at equilibrium 
and those that are not. This is accomplished by computing the critical times [4] for 
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each quasibound state. The rate constants k{ are required to calculate these critical 
times, however as the criterion for rejecting a particular critical time is approximate, 
we also approximate the rate constants by using hard sphere rate constants for the 
calculation of the critical times. These calculations are carried out for several tem- 
peratures and hard sphere parameters, and the net result is that quasibound states 
with resonance widths greater than about 0.001 cm -1 can be considered to be at 
equilibrium for the temperatures considered in the present study. Of the resonances 
characterized in Ref. [11], 38 can be considered to be at equilibrium (we use the 
resonance parameters calculated using the most accurate potential which includes 
estimates of nonadiabatic effects). In previous work [4,7,8], the list of quasibound 
states considered to be in equilibrium also included the v = 14, J = 4 state, how- 
ever more accurate calculations [ll] indicate that this state should be excluded. 
An additional restriction used previously [4,7,8] was to neglect quasibound states 
which have small equilibrium constants by virtue of their relatively high energies 
- this pared the list down to 6 quasibound states. However, because of the high 
temperatures considered here, we do not make this additional restriction. 

The three-body recombination rate constant will be evaluated using the qua- 
siclassical trajectory method. The third body M is taken to be thermal H 2 , that is 
the rate constant k{ n of Eq. (40) is given by 

k L = E (41) 

m,m' 

where m, m' are indicies labeling the bound states of H 2 , P m is the probability of 
being in state m at temperature T, and is a state-to-state rate constant. 

The probability P m is given by 


Pm — P {Jm 5 ^ m ) 

m* 


where 

P = gj(2J + l)exp(—e/ kT) 


k is Boltzmann’s constant, and 


9J 



J odd; 
J even. 


(42) 


(43) 


(44) 
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The sum in Eq. (41) is accomplished by randomly selecting the initial state m 
with probability P m as described in Ref. [28]. The states included in Eq. (41) are 
characterized using the WKB method and include only those having energies c m 
less than the dissociation energy. In particular, the energies e m are determined by 
searching for energies where the vibrational action J v given by 


is equal to 


Jv = ~ V HH {r ) - (J r 2 /2/Ztf 2 r 2 )]} 1/2 dr 

J v = {v + 1/2 )h, 


(45) 


(46) 


with v an integer. Above, nu 2 ls the reduced mass of H 2 , and J r is the classical 
rotational angular momentum and is given as 


J r = {J + l/2)h, (47) 

J the quantum mechanical rotational angular momentum quantum number. 

The sum in Eq. (40) will be evaluated in a similar manner to that used with 
Eq. (41). To do this, we write the equilibrium constant in the form 

K, = JV„(r)A(r), (48) 

where the probability P t is given by 


P, = P(J i ,t i ,T)/'£p(J,,< i .,T), (49) 

i and i 1 quasibound state indicies, and the factor N eq is given by [4] 

V 

mu the mass of an hydrogen atom. Then the recombination rate constant is given 
by 

k 3B = N eq (T) J2 Pi{T)Pm(T) Y, krm-nm>(T). (5l) 

i,m n,m # 

Each trajectory is initiated by specifying a temperature T, a maximum impact 
parameter b max , an initial separation d and by the choice of 12 random variables 
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(Ci > ^ 12 ) uniformly distributed with 0 < ft < 1 and one random variable from a 

modified normal distribution. The first random variable ft selects the initial state 
m using the probabilities P m . This determines the internal energy and classical 
rotational angular momentum J r for the bound if 2 molecule. The initial state of 
the quasibound #2 is determined by ft from the probabilities P, - this determines 
the classical rotational angular momentum. The internal energy is fixed using a 
random variable selected from a modified normal distribution as discussed in detail 
below. Next the initial relative translational energy E re i is set by finding the root 
of 

ft = 1 - exp(-E rt i/KT)( 1 4- Erei/nT), (52) 

which produces a Boltzmann distribution for E re i [28]. The initial impact parameter 
6 is then found from 

b = b max $J . (53) 

From the masses, d, b, and E re i, the initial cartesian coordinates and velocities of 
the centers of mass of the two H 2 molecules are fixed. The centers of mass are 
placed on the x axis, separated by d, and the initial velocities lie in the xz plane. 

The initial coordinates and velocities of each diatom are determined by the ran- 
dom variables ft-ft or ft-ft 2 - First the initial bond length and its time derivative are 
determined as in Ref. [29] by numerically integrating the vibrational motion over the 
fraction ft of the vibrational period. Next the orientation of the rotational angular 
momentum vector uJ is determined by the spherical polar angles arccos[2(ft — 1/2)] 
and 27rft. The magnitude of u> is given by 

w = J r /MH 2 r 2 . (54) 

Finally the rotational phase is specified by ft. With this information the cartesian 
coordinates and velocities of the atoms with respect to the center of mass of the 
diatom are fixed. The same procedure is then used with ft-ft 2 for the second 
molecule. This information, combined with the position and velocities of the centers 
of mass, completely specifies the initial conditions of the four atoms. 

We next discuss in detail our selection of the initial internal energies of the 
quasibound states. For the evaluation of the probabilities P,, we will use the en- 
ergies from the most accurate potential of Ref. [11], which includes estimates of 
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nonadiabatic effects. The use of the accurate energies presents a problem for some 
of the quasibound states, because sometimes the resonance energies lie above the 
centrifugal barrier, thus these states are not classically bound. In the previous 
studies, this problem was avoided by modifying the diatom potential[7] or by arbi- 
trarily decreasing the energy(8] so that these states became classically bound. Here 
another approach is taken. We interpret the quasibound states as states with un- 
certain energies. In particular, rather than ascribing the energy e* to all quasibound 
states labeled by t , the energies are chosen randomly from a probability distribution 
peaked about e*. The distribution used is a modified normal distribution given by 

f 0, e > Vj. - 6 ; 

P(e)de = | 0, c<<5; (55) 

( iViexp[-(^jjr i ) 2 ], otherwise, 

where Vj is the height of the centrifugal maximum for angular momentum J, 6 is 
a small number taken to be 4 /xEh, JV t - is a normalization constant, and r, is the 
resonance width for quasibound state i. The energies generated by this distribution 
are calculated using the rejection method using the unmodified normal distribution 
as the comparison function [30]. 

All uniformly distributed random numbers were generated by shuffling a multi- 
plicative congruential quasirandom number sequence. The shuffling was performed 
using algorithm B of Ref. [31] using the parameter k = 99, and the parameters 
used for the multiplicative congruential generator were modulus 2 64 and multiplier 
6364136223846793005 [31]. The reasons for using this random number generator 
are twofold. First of all it is machine independent so that calculations performed 
on different computers can easily use the same random number sequences, which is 
useful for comparing results. Secondly, and more importantly, it is a much higher 
quality random number generator than is typically available from operating system 
libraries. In the present application where we are performing monte carlo inte- 
gration of high dimensionality, it is probably important that high quality random 
numbers are used. 

The final state analysis is the same as used previously[7], namely diatoms with 
energies less than the dissociation energy are considered bound and all others are 
not. We will use stratified sampling[29] for calculating the final rate constants, i.e. 
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different trajectories for a given T will have different values of b max . We use two 
values of b ma x , called b\ and 62* Then if N% is the total number of trajectories for a 
given temperature and Nf ouni is the number of trajectories which resulted in two 
bound H 2 molecules when 1 < b < 6* (60 is zero), then the recombination rate 
constant is given by the expression 

k 3B = iMD L 8 * r - >’UW+ , ' i /N i , (56) 

*HH 3 +H 3 

Hh 3 +h 3 the reduced mass for H 2 + #2 collisions. The one sigma error associated 
with this quantity is 

(57) 

The actual integration of the trajectories will be now discussed. The integration 
coordinates used are a generalization of those used in Ref. [32], namely we initially 
require that the center of mass of the system be stationary at the origin of our 
cartesian coordinate system - then the motion of the four atoms with respect to 
the center of mass is determined by 9 cartesian coordinates and their conjugate 
momenta. We use as the 9 cartesian coordinates the coordinates of the first three 
atoms. 

In general if there are N atoms, and the center of mass is stationary at the 
origin, we can use as independent coordinates the cartesian coordinates of the first 
N — l atoms. If these coordinates are denoted by the column vector q with elements 
{ x i,yu zi,x 2 , ...,zn-i) for atoms l,2,...,iV — 1, then the kinetic energy E k i n can 
be written 

E^kin ~ 2^ ^dq, (58) 

where • denotes time derivative, T transpose, and the matrix M has elements defined 
by 

( m k + m 2 k /m N , i = j = 3 (k -1) +1,1 =1,2, or 3; 

Mij = < mkmp/mN, i = Z(k — 1) + l,j = 3(p — 1) + /,/ =1,2, or 3; (59) 

otherwise. r -~' " 

The Hamiltonian is then given by [33] 

H=^p T M~ 1 p + V, (60) 
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where V is the potential energy and the conjugate momenta p are given by 


p = Mq. (61) 

Hamilton’s equations of motion are then 


dqi/dt M l )ji 

(62) 

dpi/dt = —dV/dqi. 

(63) 


We have considered two numerical methods for integrating Eqs. (62) and (63). 
Originally we used the fourth order Adams-Bashforth-Moulton predictor-corrector 
algorithm, with fixed time steps[34], which is a common integrator used for classical 
trajectories [28]. However as pointed out previously [7], it is necessary to take very 
small time steps to obtain accurate results because at certain times the atoms can 
have very high velocities. This happens when a diatom is near its equilibrium 
separation when its kinetic energy can be close to the bond energy of H 2 . However 
most of the time, the atoms will be moving much slower, thus much computational 
effort is wasted using a fixed stepsize integrator. For this reason we seek an efficient 
variable stepsize integrator. 

The results of tests of several different integrators for classical trajectories was 
reported in Ref. [28], and the one we have chosen to develop is the Burlisch-Stoer 
method [34,35]. This method is recommended in Refs. [28,34] over several other 
methods. An important advantage which will be exploited here is that variable time 
steps can be used and programmed in a manner which does not have a detrimental 
effect on the efficient use of the vector capabilities of modern supercomputers. This 
method and our modifications to it are now described in detail. 

In the Burlish-Stoer method, in order to integrate the differential equations 
over some time increment H , one performs several calculations using a low order 
method ( the modified midpoint rule ) which differ in the time steps used. The 
time steps are from the special pattern h = (if/2, if/4, if/6, if/8, H/12 , ...), and the 
results are extrapolated to zero stepsize using a rational function. No information 
from previous steps is required for the current step, thus the method is self starting, 
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unlike the Adams-Bashforth-Moulton method, and also the time increment can be 
changed at will. It should be noted that the time interval H is typically much larger 
than the stepsizes used by other integrators, thus this method is not a good one if 
the solution of the differential equation is required on a fine time grid. 

In the implementation of integrating a single set of coupled differential equa- 
tions, it is typical to extrapolate the results after each modified midpoint calculation, 
and continue decreasing the time step until the extrapolation converges. Then the 
time interval H is adjusted so that the expected number of integrations required for 
convergence at the next time step will be some predetermined optimum number [34]. 
However in the present application, we will take advantage of vector processing by 
simultaneously integrating a number of independent trajectories, thus it will not be 
efficient to keep decreasing the time step until the trajectory which is most diffi- 
cult to integrate is converged. Thus we modify the above scheme by integrating all 
trajectories using a predetermined number of stepsizes. Each trajectory will have 
its own time increment. After the integrations over a time increment are carried 
out, the errors in the extrapolations from the various stepsizes are checked. If the 
errors are small enough, the time interval H for that trajectory is increased and 
the clock is advanced. If the errors are too large, the results from the time incre- 
ment are thrown away, H is reduced and the clock is not advanced. Occasionally 
the trajectories are checked for completion, and if required, finished trajectories are 
removed and new ones are started. This scheme is vectorizable except where the 
time increment H is changed and end tests are performed, but these steps usually 
amount to an insignificant portion of the total resources required in a calculation, 
thus this algorithm runs efficiently on vector pipeline supercomputers. 

We now discuss some details of the algorithm, in particular how the errors are 
estimated and the time interval is changed. Let the number of different time steps 
per time interval H be N TS and define the quantity ^ to be the j’th integration 
function (one of the 6 N — 6 coordinates and conjugate momenta) resulting from 
the *’th extrapolation using the first i time steps. We estimate the errors of the 
extrapolations to be 

N TS -l. (64) 
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To remove the units in this error and to smooth out the variations from trajectory 
to trajectory, this error is scaled by the rms initial value of either all coordinates or 
all momenta for this trajectory - this produces a relative error called Afj^. Finally 
the maximum error over all j is called A^l and this is used to characterize the error 
obtained using t time steps. This error is compared to an error tolerance parameter 
e, and the largest value of i for which A^ 1 ^ is greater than e is determined - this 
value of i is called i c — 1. If all A^ are greater than or equal to £, then the 
error is unacceptably large and the time increment H is reduced by a factor of 
l/< T e ? /2 , where Nf Uv is the ratio of the time increment H to the time step h for 
integration i. If i c is less than N TS , then the time increment is increased by the 
factor {{N^rs I + iu]/(l 4- w), where w is a fixed damping parameter which 
we usually take equal to 1/2. 

In practice, this time increment changing algorithm works quite well - see 
Fig. (4) where the relative error versus CPU time is compared for the current 
method and the Adams-Bashforth-Moulton method. For the present application 
the Burlisch-Stoer method is always more efficient then the Adams method, and 
its superiority increases rapidly as the desired relative error decreases. Thus there 
is less motivation to scrimp using the present method, for a small increase in CPU 
time can produce a large increase in accuracy. The data for this figure is based on 
calculations with N TS equal to 5, w equal to 1/2, which experiments show to be 
about optimum for the present application. In general larger values of N TS allow 
larger time increments but if the time increment becomes too large, the terrain 
covered in a single increment will be sufficiently varied so that it is not possible to 
make significant adjustments to the time increment. This tends to favor smaller 
values of N TS for the present application. Initial conditions for the trajectories 
were selected as described above for a temperature of 5000 K, b m ax equal to 8 a 0 , 
and an initial separation d equal to 15 a 0 . The potential V jj° n was use for these 
calculations. Ten trajectories were integrated forward for a time of 8000 a.u., then 
back integrated for a time 8000 a.u. The times reported as the abscissa are the 
average CPU seconds for each trajectory using a single CPU on the Ames ACF 
Cray-XMP/48. Since these were obtained running only ten trajectories and did 
not use the accurate potential, they are not indicative of the times required for 
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our production calculations. We estimate that our production calculations require 
about a factor of five less time under these conditions. The relative error reported 
is the rms relative difference between all of the initial position and momentum 
variables and their back integrated values. The value of 8000 a.u. is a typical 
integration time for the conditions used here. The initial value of H was take to be 
10 a.u. for all calculations of Fig. (4). As an example of the changes made to H, 
when e was lxlO -5 , H ranged from 7.6 to 87.7 a.u., an average of 195 successful 
and 23 unsuccessful time increment steps were taken per trajectory and the rms 
relative error was 6.2 xlO -4 . Another advantage of a variable stepsize integrator 
is that the errors for each trajectory are about the same - when the rms relative 
error was 6.2 xlO -4 , the rms relative errors for individual trajectories varied by 
a factor of 53, while when using the Adams-Bashforth-Moulton method with a 
stepsize giving a similar rms relative error of 1.9 xlO -4 , the rms relative errors for 
individual trajectories varied by a factor of 3200. 

Although the increment changing algorithm works well, it suffers from sev- 
eral potential drawbacks. First consider the determination of the error A(j*\ 
As presently described, if one were to change the initial separation between the 
molecules, but nothing else, the integrator would allow larger errors because the 
scaling factor, which is based on the rms initial conditions, would be increased. An- 
other manifestation of this is that trajectories with larger initial impact parameters 
would be allowed to have larger errors. An additional question is what the relative 
weighting of the errors in the coordinates and the momenta should be for optimum 
performance. Turning to the time increment changing algorithm, it always tries to 
increase the time increment until it is too big and a step is rejected. It would be 
better to have some method of slowly decreasing the time increment as the errors 
get larger, however this would increase the complexity of the algorithm. Once it 
has been decided to increase the time increment, the factor we use to increase it is 
somewhat arbitrary, and situations may arise where increasing the time increment 
by Ntffs /Nfc+i would lead to a excessively large increment. The inclusion of the 
damping parameter w is an attempt to forestall this possibility of increasing the 
time increment too quickly and then having to decrease it right away. The time in- 
crement decrease factor is also arbitrary, however it seems to satisfy two important 
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properties, namely that it is unlikely that two decreases occur consecutively and 
that the increase algorithm can overcome the rapid decrease in a few applications, 
if the error is small enough. 

An additional aspect of the present calculations concerns the question of the 
number of trajectories to be integrated simultaneously. The larger the number, the 
more efficient the calculations will be since we use vector pipeline supercomputers 
(various calculations were performed the the Ames ACF Cyber 205, Cray-XMP/48, 
or the NAS facility Cray-2). On the other hand, larger numbers require more com- 
putational resources. For most production calculations, we integrate a maximum 
of 500 trajectories simultaneously. This number is large enough so that significant 
additional speedups are not possible by increasing it further. 

IV. Results and Discussion 

The main results of this paper are given in Table IX. The rate constants for each 
temperature are based on the results of 4000 trajectories which produces statistical 
one sigma errors of about 5 %. All calculations use an initial separation d equal 
to 15 a 0 and were terminated when any atom-atom distance exceeded 27 a 0 . The 
maximum impact parameters (61,62) varied from (6.4,8) a 0 at 5000 K to (8.4,12) 
a 0 at 100 K. Most calculations used the integration parameters N TS = 5, w = 1/2, 
e = 1 x 10 -5 , and a initial time increment of 10 a.u. The remaining calculations 
used an earlier version of the error control algorithm and should produce results of 
commensurate accuracy. 

At low temperatures, no exchange processes occurred and so it was probably 
not necessary to force the potential to be symmetric with respect to exchanging 
atoms, but at higher temperatures this symmetry is more important. For 2000 K, 
0.13 % of the trajectories ended up in other than the initial arrangement AB + CD, 
at 3000 K the number was 1.0 % and at 5000 K, 3.7 %. Also it is not known 
how many trajectories crossed only temporarily into other arrangements, thus it 
is probably good that our potential is symmetric. The 144 kcal/mol barrier to 
exchange may seem too large to allow any exchange reactions, but the quasibound 
states have at least 109 kcal/mol of energy making the effective barrier no more 
than 35 kcal/mol. Thus the exchange rates we observe are not unreasonable. 

In addition to Table IX, our results are also shown in Fig. 5 along with several 
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experimental measurements of the recombination rate constantsfl, 36-41], We show 
the results of the measurements of one laboratory for a range of low temperatures 
where the measurements of various groups agree fairly well, and the results of 
several groups at high temperature to show the spread obtained there. At 100 
K, the theoretical results are too small by a factor of about 3, and at the higher 
temperatures, the theoretical results are about half of the recommended values. 
The temperature dependence of the present results are similar to estimates of the 
experimental temperature dependence, however they tend to fall off too quickly 
at low temperature ( T < 200 K). For high temperature, the theoretical results 
also appear to fall off too quickly, however the situation is less clear because of the 
uncertainty in the experimental measurements. In analyzing the experiments, it is 
common to assume that k 3B = A/T m [36-38,40], and some in experiments, there 
is insufficient information to accurately determine the exponent m so a value of 
one is assumed [36-38]. Thus the experimental temperature dependence is not well 
determined over the range 3000-5000 K. 

Because of the relatively poor agreement with experiment, it is interesting to 
compare to the results of previous calculations. The two temperatures in common 
are 100 and 300 K where we obtain 1.84 xlO 15 cm 6 mol ~ 2 sec -1 and 1.52 xlO 15 , 
Whitlock et al.[ 7] obtain 6.6 xlO 15 and 4.2 xlO 15 using their DHR potential, and 
experiment gives 6. 2±0.7 x 10 15 (this is estimated from the graph of Ref. [41]) 
and 3.0=0. 2 x 10 15 ( at 298 K [41]). The results of Whitlock et al. [7] are in 
better agreement with experiment than our presumably more accurate calculations. 
A detailed analysis shows that the differences between the two calculations are 
primarily due to a combination of the differences in the interaction potential and 
the differences in the diatomic potential. If we repeat our calculations using the 
DHR potential of Ref. [7] instead of the ab initio interaction potential of Sect. II 
(we use our accurate diatomic potential), we obtain 4.09 xlO 15 and 2.73 xlO 15 . 
The 300 K results of this calculation are in even better agreement with experiment 
and further calculations using the DHR potential at higher temperatures continue 
this trend [see Fig. (6)]. 

We may then ask if the DHR potential is more accurate than the ab initio 
potential since it produces results which are in closer agreement with experiment. 
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If this were so, it would be quite surprising for the DHR potential is very simple, 
comprising of the sum of pairwise interactions. It is not easy to pin point the features 
of the potentials primarily responsible for the differences in the dynamics, however 
it is possible to single out two differences. First of all consider Fig. (2) where the 
dependence of the interaction potential on bond length is shown. The ab initio 
calculations clearly indicate that for bond lengths near equilibrium, the potential 
increases as the bond length increases and as the bond length becomes larger, this 
trend reverses. In contrast, for the geometry of this figure, any pairwise potential, 
including the DHR potential, will cause the potential to monotonically decrease as 
the bond length increases. This difference in force should have a profound effect 
when both molecules are near their equilibrium displacement [42]. Next consider the 
anisotropies of the two potentials. These are also very different when both molecules 
are near their equilibrium separation - this is shown in Fig. (3). For R less than 5 a 0 , 
the DHR potential gives a leading anisotropy which is an order of magnitude larger 
than the ab initio potential. The spherical average of the DHR potential agrees well 
with the ab initio potential. The differences in the anisotropies should have a big 
effect on rotational energy transfer. It is beyond the realm of possibility that the 
ab initio calculations have errors so large that the DHR potential is the accurate 
potential for the features just discussed, however these features are not necessary 
the primary progenitors of the differences in the dynamics calculations. However, 
it is most likely that the ab initio potential is much more globally accurate than the 
DHR potential and that the good agreement with experiment obtained using the 
DHR potential is due to fortuitous cancellations of error. 

There are several possible reasons why the present calculations underestimate 
experiment by such a large factor. First of all there are errors in the present results 
because of the use of an approximate potential energy function. However prelim- 
inary results obtained on different potentials based on the ab initio calculations 
were quite similar to the present results and so it seems unlikely that reasonable 
variations in the potential alone could make up for the difference. Another fault 
of our calculations is that they use classical dynamics. Classical dynamics has not 
been tested for a system such as this one, however one usually expects reasonable 
accuracy for calculations involving highly excited states, such as those involved in 


29 



the present study. On the other hand, classical dynamics is unreliable for describ- 
ing processes for which tunneling plays an important role. Thus the present use 
of classical dynamics must be treated with some suspicion because of the presence 
of hydrogen atoms which can undergo large amounts of tunneling. In addition, 
tunneling plays an important role in the formation and decay of the quasibound 
states. 

Furthermore, several processes are omitted from the present calculations, 
namely the possibility of direct association from the continuum as well as contribu- 
tions to the recombination rate from other processes involving quasibound states, 
e.g. the chaperon mechanism: 


H + H 2 ~H' 3 , 

(65) 

H~Hz^ 2H 2 . 

(66) 


Estimates from Ref. (4) indicate that direct association from the nonresonant con- 
tinuum relative to the association from the quasibound states will become more 
important as the temperature increases, however if this was the main cause for the 
underestimate, then one would expect that the present results would underestimate 
experiment by increasing amounts as the temperature increases. This contrasts to 
the approximately constant underestimate of a factor of two over a fairly broad 
temperature range, thus it seems unlikely that the neglect of direct association is 
the primary fault of the present calculations. As far as the chaperon mechanism is 
concerned, it seems to be important only at lower energies than we consider here 
[5-7] . However, this mechanism may be responsible for the underestimate at low 
temperatures [see Figs. (5)-(6)|. It should be noted that a converged quantum 
mechanical dynamics calculation would properly take into account the role of the 
continuum and the chaperon mechanism 

Finally one can consider the possibility that the resonance complex mechanism 
is not the proper route for calculating the recombination rate constants - this is 
because the bottleneck for the association may occur elsewhere when nonequilibrium 
and redissociation processes are considered. In particular Ashton et al. [43] conclude 
from their model calculations that the primary bottleneck occurs in the relaxation 
of bound vibrational states close to the dissociation limit. Thus the results of 
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the resonance complex mechanism should be an upper bound to the true rate of 
recombination which then means that the underestimate of the present calculations 
is even greater than indicated in Fig. (5). 

There exists one more possibility for explaining the underestimate, namely 
that the experimental results are systematically high. However this seems unlikely 
due to the degree of agreement between various experimental workers and we will 
not consider this as an serious proposition until the approximations in the current 
theoretical work have been more thoroughly tested. 

V. Conclusions 

We have carried out calculations of the recombination rate constants for the 
process H + H + H 2 — » ► 2 #2 using a quasiclassical trajectory implementation of the 
energy transfer mechanism of the resonance complex theory. A new global poten- 
tial energy function for H 4 is constructed which should be accurate for geometries 
important to the dynamics calculations. In spite of the quality of these inputs into 
the calculation, the agreement with experimental results is not very good with the 
theoretical results about a factor of 2 too small over the temperature range 300 - 
5000 K. Of the are several possible causes for this underestimate, the most likely 
reason is the importance of quantum effects. This and other limitations must be 
studied in order to predict with confidence three body recombination rate constants 
for hydrogen containing compounds. 

VI. Acknowlegdements 

Financial support for this work is through NASA grant NCC 2-443. Calcula- 
tions on the Cray-2 were made possible by a grant of computer time from the NAS 
facility. In addition, we wish to thank C. W. Bauschlicher for assistance with the 
electronic structure calculations. 



References: 

1. N. Cohen and K. R. Westberg, J. Phys. Chem. Ref. Data 12, 531 (1983). 

2. D. L. Bunker, J. Chem. Phys. 32, 1001 (1960). 

3. J. C. Keck, J. Chem. Phys. 32, 1035 (1960). 

4. R. E. Roberts, R. B. Bernstein, and C. F. Curtiss, J. Chem. Phys. 50, 5163 
(1969). 

5. V. H. Shui and J. P. Appleton, J. Chem. Phys. 55, 3126 (1971). 

6. R. T Pack, R. L. Snow, and W. D. Smith, J. Chem. Phys. 56, 926 (1972). 

7. P. A. Whitlock, J. T. Muckerman, and R. E. Roberts, Chem. Phys. Lett. 16, 
460 (1972); J. Chem. Phys. 60, 3658 (1974). 

8. A. E. Orel, J. Chem. Phys. 8T, 314 (1987). 

9. R. J. LeRoy- and R. B. Bernstein, J. Chem. Phys. 54, 5114 (1971). 

10. C. Schwartz and R. J. LeRoy, J. Mol. Spect. 121, 420 (1987); R. J. LeRoy 
and C. Schwartz, University of Waterloo Chemical Physics Research Report CP-301 
(1986). 

11. D. W. Schwenke, Theor. Chim. Acta, submitted. 

12. W. Kolos, Int. J. Quan. Chem. 1, 169 (1967). 

13. W. J. Deal and R. H. Young, Mol. Phys. 19, 427 (1970). 

14. W. Kolos, K. Szalewicz, and H. J. Monkhorst, J. Chem. Phys. 84, 3278 (1986). 

15. MOLECULE is a vectorized gaussian integral code written by J. Almlof and 
SWEDEN is a vectorized SCF-MCSCF-direct Cl -conventional CI-CPF-MCPF pro- 
gram written by P. E. M. Siegbahn, C. W. Bauschlicher, B. Roos, P. R. Taylor, A. 
Heiberg, J. Almlof, S. R. Langhoff, and D. P. Chong. 

16. F. B. Van Duijenveldt, IBM Tech. Res. Rep. RJ947 (no. 16437) (1979). 

17. M. J. Frisch, J. A. Pople, and J. S. Binkley, J. Chem. Phys. 80, 3265 (1984). 

18. W. Meyer, P. C. Hariharan, and W. Kutzelnigg, J. Chem. Phys. 73, 1880 
(1980). 

19. G. Gioumousis and C. F. Curtiss, J. Math. Phys. 2, 96 (1961). 

20. W. Kolos and L. Wolniewicz, J. Chem. Phys. 43, 2429 (1965). 

21. M. J. Norman, R. O. Watts, and U. Buck, J. Chem. Phys. 81, 3500 (1984). 

22. K. T. Tang, J. M. Norbeck, and P. R. Certain, J. Chem. Phys. 64, 3063 (1976). 


32 



23. F. B. Brown, D. W. Schwenke, and D. G. Truhlar, Theor. Chim. Acta 68, 23 
(1985). 

24. B. Liu, J. Chem. Phys. 58, 1924 (1973); P. Siegbahn and B. Liu, J. Chem. 
Phys. 68, 2457 (1978); D. G. Truhlar and C. J. Horowitz, J. Chem. Phys. 68, 2466 
(1978); 71, 1514 (1979). 

25. P. G. Burton and U. E. Senff, J. Chem. Phys. 76, 6073 (1982). 

26. M. Rubinstein and I. Shavitt, J. Chem. Phys. 51, 2014 (1969). 

27. D. M. Silver and R. M. Stevens, J. Chem. Phys. 59, 3378 (1973). 

28. D. G. Truhlar and J. T. Muckerman, in Atom-Molecule Collision Theory, ed. 
by R. B. Bernstein (Plenum Press, New York, 1979), pg. 505. 

29. N. C. Blais and D. G. Truhlar, J. Chem. Phys. 65, 5335 (1976). 

30. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical 

Recipes (Cambridge University Press, Cambridge, 1986), Chapter 7. 

31. D. E. Knuth, Seminumerical Algorithms, Vol. 2 of The Art of Computer Pro- 
gramming, (Addison- Wesley, Reading Mass., 1981), 2nd. edition, Chapter 3. 

32. D. L. Bunker and N. C. Blais, J. Chem. Phys. 41, 2377 (1964). 

33. H. Goldstein, Classical Mechanics (Addison- Wesley, Reading Mass., 1980), 2nd 
e<3(cfition, pg. 344. 

34. Reference 30, Chapter 15. 

35. R. Bulirsch and J. Stoer, Numer. Math. 8, 1 (1966). 

36. J. P. Rink, J. Chem. Phys. 36, 262 (1962). 

37. R. W. Patch, J. Chem. Phys. 36, 1919 (1962). 

38. E. A. Sutton, J. Chem. Phys. 36, 2923 (1962). 

39. T. A. Jacobs, R. R. Giedt, and N. Cohen, J. Chem. Phys. 47, 54 (1967). 

40. I. R. Hurle, A. Jones, J. L. J. Rosenfeld, Proc. Roy. Soc. A 310, 253 (1969). 

41. D. O. Ham, D. W. Trainor, and F. Kaufman, J. Chem. Phys. 53, 4395 (1970). 

42. J. W. Duff and D. G. Truhlar, J. Chem. Phys. 63, 4418 (1975). 

43. T. Ashton, D. L. S. McElwain, and H. O. Pritchard, Can. J. Chem. 51, 237 
(1973). 


33 


Table I. Parameters for the analytic representation given in Eqs. (l)-(4) of the 
1 H 2 X j potential curve. All parameters are in atomic units and are valid only 

for the isotope of H with mass 1 amu. 


r 0 

3.5284882 


-1.272679684173909 

Cl 

6.499027 

O 7 

0.5630423099212294 

Q 

124.3991 

08 

-0.1879397372273814 

r a 

'-'10 

3285.828 

O 9 

0.04719891893374140 

D 

0.160979391 

010 

-0.008851622656489644 

s 

1 

Oil 

0.001224998776243630 

CL 0 

0.03537359271649620 

Oi2 

-1.227820520228028(-4) 6 

a x 

2.013977588700072 

Ol3 

8.638783190083473(— 6) 

02 

-2.827452449964767 

Ol4 

— 4.036967926499151(— 7) 

®3 

2.713257715593500 

®15 

1.123286608335365(— 8) 

O 4 

-2.792039234205731 

Ol6 

— 1.406619156782167(— 10) 

05 

2.166542078766724 




a From Ref. [13]. 

b The number in parenthesis is the power of ten multiplying the number. 


Table II. Computed energies in Eh for separated H? + H2 . 


r 1 

T2 

gCAS 

EMRci 

1.000 


- 2.2501723 

- 2.2951620 

1.401 

■SSH 

- 2.3040214 

- 2.3456759 

3.000 

mMSSm 

- 2.2002567 

- 2.2286053 


1.401 

- 2.1548434 

- 2.1762882 

20.00 

1.401 

- 2.1519998 

- 2.1729908 

1.000 

1.000 

- 2.1963232 

- 2.2446483 

20.00 

20.00 

- 1.9999781 

- 1.9999782 


Table III. Computed energies in Eh for linear geometries. 


r i 

T2 

R 

gCAS 

gMRCI 


1.000 

3.0 

- 2.1489968 

- 2.2021605 

1.401 

1.000 

3.0 

- 2.1895079 

- 2.2411766 

1.401 

1.401 

3.0 

- 2.2277344 

- 2.2782835 


1.000 

3.5 

- 2.1769510 

- 2.2280992 


1.000 

3.5 

- 2.2242814 

- 2.2732647 

1.401 

1.401 

3.5 

- 2.2708059 

- 2.3180081 


1.000 

4.0 

- 2.1885540 

- 2.2384542 

1.401 

1.000 

4.0 

- 2.2392485 

- 2.2864936 

1.401 

1.401 

4.0 

- 2.2894001 

- 2.3342640 

1.000 

1.000 

4.5 

- 2.1932777 

- 2.2424759 

1.000 

1.401 

4.5 

- 2.2456559 

— 2,2918954 


1.000 

4.5 

— 2.2456556 

- 2.2918956 

1.401 

1.401 

4.5 

- 2.2976739 

- 2.3411200 

3.000 

1.401 

4.5 

- 2.1810055 

- 2.2164670 

5.000 

1.401 

4.5 

- 2.0795631 

- 2.1188059 

1.000 

1.000 

5.0 

- 2.1951517 

- 2.2439630 

1.401 

1.000 

5.0 

- 2.2483393 

- 2.2940230 


1.401 

5.0 

- 2.3013097 

- 2.3439591 






Table IV. Computed energies in Eh for T shaped geometries. 


r l 

T2 

R 

£CAS 

jgMRCI 

1.000 

1.000 

3.0 

-2.1581148 

-2.2109995 

1.000 

1.401 

3.0 

-2.2026122 

-2.2541250 

1.401 

1.000 

3.0 

-2.2083495 

-2.2584915 

1.401 

1.401 

3.0 

-2.2534653 

-2.3025120 

1.000 

1.000 

3.5 

-2.1806088 

-2.2315200 

1.000 

1.401 

3.5 

-2.2294472 

-2.2781860 

1.401 

1.000 

3.5 

-2.2322244 

-2.2802336 

1.401 

1.401 

3.5 

-2.2811433 

-2.3271734 

3.000 

1.401 

3.5 

-2.1774572 

-2.2105028 

5.000 

1.401 

3.5 

-2.1459055 

-2.1699190 

1.000 

1.000 

4.0 

-2.1900624 

-2.2398409 

1.000 

1.401 

4.0 

-2.2414512 

-2.2885466 

1.401 

1.000 

4.0 

-2.2427205 

-2.2894577 

1.401 

1.401 

4.0 

-2.2940422 

-2.3382245 

1.401 

3.000 

4.0 

-2.1715427 

-2.2091528 

3.000 

1.401 

4.0 

-2.1892140 

-2.2205532 

1.401 

5.000 

4.0 

-2.0467742 

-2.0870823 

5.000 

1.401 

4.0 

-2.1500841 

-2.1731692 

1.000 

1.000 

4.5 

-2.1939053 

-2.2430450 

1.000 

1.401 

4.5 

-2.2466281 

-2.2927943 

1.401 

1.000 

4.5 

-2.2471777 

-2.2931716 

1.401 

1.401 

4.5 

-2.2998260 

-2.3429269 

1.000 

1.000 

5.0 

-2.1954200 

-2.2442018 

1.000 

1.401 

5.0 

-2.2487844 

-2.2944316 

1.401 

1.000 

5.0 

-2.2490111 

-2.2945755 

1.401 

1.401 

5.0 

- 2.3023270 

- 2.3448057 


Table V. Computed energies in Eh for parallel shaped geometries. 


r 1 

r 2 

R 

jgCAS 

gMRCI 

1.000 

1.000 

3.00 

- 2.1630376 

- 2.2141120 

1.401 

1.000 

3.00 

- 2.2125055 

- 2.2605994 

1.401 

1.401 

3.00 

- 2.2614119 

- 2.3065803 

1.401 

1.401 

3.00 

- 2.2468856 

- 2.3062874 

1.000 

1.000 

3.50 

- 2.1827285 

- 2.2326254 

1.401 

1.000 

3.50 

- 2.2340995 

- 2.2809217 

1.401 

1.401 

3.50 

- 2.2851299 

- 2.3289184 

1.000 

1.401 

3.75 

- 2.2397741 

- 2.2861676 

1.401 

1.401 

3.75 

- 2.2915579 

- 2.3348676 

3.000 

1.401 

3.75 

- 2.1853072 

- 2.2155281 

5.000 

1.401 

3.75 

- 2.1486507 

- 2.1716260 

1.000 

1.000 

4.00 

- 2.1909127 

- 2.2401360 

1.401 

1.000 

4.00 

- 2.2434838 

- 2.2895480 

1.401 

1.401 

4.00 

- 2.2958436 

- 2.3387790 

1.000 

1.000 

4.50 

- 2.1942165 

- 2.2430549 

1.401 

1.000 

4.50 

- 2.2474464 

- 2.2930622 

1.401 

1.401 

4.50 

- 2.3005503 

- 2.3429644 

1.000 

1.000 

5.00 

- 2.1955172 

- 2.2441340 

1.401 

1.000 

5.00 

- 2.2490789 

- 2.2944289 

1.401 

1.401 

5.00 

- 2.3025690 

- 2.3446655 


38 


Table VI. Computed energies in Eh for cross shaped geometries. . ■ ■ 


r 1 

T 2 

R 

gCAS 

gMRCI 

1.401 

1.401 

2.00 

- 2.1169936 

- 2.1692110 

mm 


3.00 

- 2.1631831 

- 2.2146856 

EH 


3.00 

- 2.2129581 

- 2.2615710 

1.401 

1.401 

3.00 

- 2.2624792 

- 2.3082718 

1.000 


3.50 

- 2.1827569 

- 2.2328906 

1.401 


3.50 

- 2.2342374 

- 2.2813539 



3.50 

- 2.2854980 

- 2.3296479 

1.000 

1.401 

3.75 

- 2.2398544 

- 2.2864620 

1.401 

1.401 

3.75 

- 2.2917836 

- 2.3353583 

3.000 

1.401 

3.75 

- 2.1861503 

- 2.2165969 



3.75 

- 2.1491633 

- 2.1719592 

1.000 

1.000 

4.00 

- 2.1909183 

- 2.2402610 

1.401 

1.000 

4.00 

- 2.2435329 

- 2.2897484 

1.401 


4.00 

- 2.2959868 

- 2.3391111 


1.000 

4.50 

- 2.1942185 

- 2.2431111 

1.401 

1.000 

4.50 

- 2.2474683 

- 2.2931541 

1.401 

1.401 

4.50 

- 2.3006153 

- 2.3431180 

1.000 

1.000 

5.00 

- 2.1955189 

- 2.2441593 

1.401 

1.000 

5.00 

- 2.2490912 

- 2.2944721 

1.401 

1.401 

5.00 

- 2.3026036 

- 2.3447394 

1.401 

1.401 

6.00 

- 2.3037829 

- 2.3455943 
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Table VII. Parameters for the analytic representation of the interaction potential. 
All parameters are in atomic units. 


OQ 

0.685147 

mg 

3 

b Q 

- 0.198019 

ag 

0.196784 

C Q 

0.144382 

H 

0.508106 

(Lq 

0.449685 

d\ R 

54 

mg 

1 

4 * 

0.4 

ag 

1.17467 

A p 

0.717277143 

^6 

1.72657 

B P 

1.00471655 

ml 

2 

c P 

1.27158465 

a\ 

0.0198867 

3 a 

7 

H 

0.0681996 

(3 

0.1 

3 

OOO 

2 

Po 

11 

ag 

9.52859 

a D 

2 


0.796373 

r D 

5.4 



0 A 

20 
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Table VIII. Parameters for the multiplicative correction factor of Eq. (17) a . All 
parameters are in atomic units. 


ijk 

c' jk 

c ooo 

Ji k 

c 200 

c '3 k 

c 220 

c 222 

0 0 0 

8.2335075(-l) 6 

9.4069522(— 2) 

7.914283l(— 3) 

— 4.7369816(— 3) 

0 1 0 

1.9927182(— 1) 

— 2.6657323(— 2) 

— 8.7591868(— 3) 

9.4419462(— 3) 

10 0 

1.9927182( — 1) 

— 1.7873925(— 1) 

— 8.7591868(— 3) 

9.4419462(— 3) 

110 

— 2.7306801(— 1) 

7.6910059(— 2) 

5.4783229(— 3) 

— 1.6362036(— 2) 

0 0 1 

— 2.1105535(— 1) 

— 1.7750047(— 2) 

4.5662762(— 4) 

5.3064904(— 4) 

0 1 1 

8.1872059(— 3) 

1.6826832(— 2) 

— 2.4577906(— 3) 

— 1.7568388(— 3) 

10 1 

8.1872059(— 3) 

3.7743965(— 2) 

— 2.4577906(— 3) 

— 1.7568388(— 3) 

1 1 1 

8.1615668(— 2) 

— 3.8437296(— 2) 

8.2469313(— 3) 

4.0949053(— 3) 


ij k Jik 

C 020 — ^OO- 


The number in parenthesis is the power of ten multiplying the number. 
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Table DC. Three body recombination rate constants. 


T(K) 

k 3B (10 14 cm 6 mol 2 sec ~ 1 ) 

100 

18.4 ± 4.5% a 

300 

15.2 ± 5.4% 

1000 

9.82 ± 5.4% 

2000 

6.81 ± 6.2% 

3000 

4.87 ± 6.6% 

5000 

3.31 ± 6.9% 


a One sigma error bars. 
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Figure Captions. 

Fig. (1): Effective Morse exponent a and r times a as a function of r for the H 2 
potential. 

Fig. (2): The H 2 + H 2 interaction potential for R = 4, ri = 1.4 a 0 , 0\ = <j> = 0, 
= 7t/ 2 as a function of r 2 for various stages of the fit. 

Fig. (3): Angular expansion coefficients v qiq2fi of the interaction potential as a 
function of R. The bond lengths are 1.449 a 0 . The DHR potential is from Ref. [7] 
and the PNO-CI and CEPA2 results are from Ref. [25]. 

Fig. (4): A comparison of the relative errors obtained using the Burlisch-Stoer 
(B.S.) algorithm and the Adams-Bashforth-Moulton (A.B.M.) algorithm as a func- 
tion of CPU time. 

Fig. (5): A comparison of the theoretical and experimental recombination rate 
constants. 

Fig. (6): A comparison of the results using the potential of Sect. II and the results 
obtained using the DHR potential of Ref. [7] with the estimate of the recombination 
rate based on experiment [l] . 
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